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ABSTRACT 

Context. Analytic methods to investigate periodic orbits in galactic potentials. 

Aims. To evaluate the quality of the approximation of periodic orbits in the logarithmic potential constructed using perturbation theory 

based on Hamiltonian normal forms. 

Methods. The solutions of the equations of motion corresponding to periodic orbits are obtained as series expansions computed by 

inverting the normalizing canonical transformation. To improve the convergence of the series a resummation based on a continued 

fraction may be performed. This method is analogous to that looking for approximate rational solutions (Prendergast method). 

Results. It is shown that with a normal form truncated at the lowest order incorporating the relevant resonance it is possible to 

construct quite accurate solutions both for normal modes and periodic orbits in general position. 

Key words, galaxies: kinematics and dynamics ~ methods: analytical 



> 

cn 
o\ 
i> 

(N 

O 
oo 
O 



X 



1. Introduction 



In his book on Dynamical Astronomy, Contopoulos (2004) en- 
courages to investigate higher-order versions of the Prendergast 
d 19821 1 method to solve non-linear differential equations. The 
original method was applied by Contopoulos & Seimenis ( 119901 
hereafter CS90) to periodic orbits in the logarithmic potential 
and consists in approximating the exact solution with rational 
trigonometric functions. Even though the trigonometric series 
used in the rational approximation are truncated at the first non- 
trivial order, in CS90 is shown that the quality of the fit to the 
exact result is quite good over a wide range of energy and el- 
lipticity. On this basis it is natural to presume that higher-order 
truncations would improve the quality of the prediction. 

However, even the simplest version of the Prendergast 
( 1198 2') method has two problematic aspects: 1) the choice of the 
dominant harmonic in the trigonometric series has to be made on 
the basis of some knowledge about the orbit type under study; 2) 
the determination of the coefficients in the series, which depend 
on the parameters of the system and on initial conditions, orig- 
inates from a non-linear algebraic system the solution of which 
must in general be performed numerically. This second aspect 
spoils the approach of much of its simplicity, even more if we 
aim at higher order truncations and consider the growth of the 
number of unknown coefficients. 

In this paper we would like to explore the link between the 
Prendergast-Contopoulos approach and the approximation of or- 
bital solutions found with a resonant normal form. The motiva- 
tion for this study stems from the idea of rooting a simplified 
version of the rational solution method into the frame of a modi- 
fied normalization algorithm in order to devise a completely an- 
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alytical approach. In fact it has recently proposed (Pucacco et al. 
[2008J to exploit a resummation technique based on continued 
fractions to speed up the convergence of series obtained in the 
framework of normal form perturbation theory. This technique 
is able to extend the quality of predictions concerning the insta- 
bility of normal modes and consequent bifurcations of families 
of boxlets (Belmonte et al. 120071) . 

In analogy with CS90 we apply this approach to investigate 
periodic orbits in the logarithmic potential (Binney & Tremaine 
1987 ). We find analytical solutions of the equations of motion for 
the normal modes and the main low-order boxlets ('loops' and 
'bananas'). By inverting the normalizing transformation of co- 
ordinates, these solutions are either in the form of standard trun- 
cated power series or in a rational form constructed by a contin- 
ued fraction truncated at the same order of the series. Knowing 
the 'normal form' approximating the system under study, the 
procedure of creating those solutions is straightforward and does 
not require any further approximation or numerics. 

We show that the analytic rational solutions obtained in this 
way offer a degree of reliability comparable, where data are 
available, to those of the semi-analytic treatment based on the 
Prendergast-Contopoulos approach. Both loops and bananas are 
quite well reconstructed in shape and dimension. We extend the 
analysis in CS90 also to check the energy conservation along 
the boxlets: it turns out that, whether for normal modes energy 
is conserved within a few percent, for loops and bananas, at this 
level of approximation, it is not easy to go below 10%. 

The plan of the paper is as follows: in Section 2 we briefly 
recall the method to construct normal forms for the logarithmic 
potential, relegating to the Appendix the explicit expressions of 
the 1:1 and 1:2 Hamiltonian and generating function. In Section 
3 we analyze the approximation of the major-axis orbit and in 
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Sections 4 and 5 we do the same respectively for the loop and 
banana families. In Section 6 we sketch the conclusions. 



2. Normal forms for the logarithmic potential 

We investigate the dynamics in the potential 



V^^log{R^ + x^ 



(1) 



For every finite values of the "core radius" R, the choice R = I 
can be done without any loss of generality. However, it is also 
of relevance the singular limit R -^ associated to a central 
density cusp (Miralda-Escude & Schwarzschild ll989b . With the 
choice R = I, the energy E may take any non-negative value. 
Otherwise, the singular limit is 'scale-free' and the dynamics are 
the same at every energy. The parameter q gives the "ellipticity" 
of the figure and ranges in the interval 



0.6 <^< 1. 



(2) 



Lower values of q' can in principle be considered but coiTespond 
to an unphysical density distribution. Values greater than unity 
are included in the treatment by reversing the role of the coordi- 
nate axes. 

Normal forms for the Hamiltonian system corresponding 
to the potential ^ are constructed with standard methods 
(Boccaletti & Pucacco 1999 Giorgilli [2002l) and have been used 
to determine the main features of its orbit structure (Belmonte et 
al. 2006, 2007 ). The starting point is the series expansion of ([TJ 
around the origin 



1 



2 2-2 
where 



q^ 



1 2 1 

s + 



2-3 



■s'- 



1 
2~4^ 



./ + . 



(3) 



(4) 



Here we briefly resume the procedure just with the purpose of 
fixing notations. After a scaling transformation 



Py — > yfqpy, y — ^yl^lq, 

the original Hamiltonian 

H(p_„py, x,y) = -(pI + pl/q) + V(s(x,y)) 



(5) 



(6) 



is subject to a canonical transformation to new variables 
Px,Py,X,Y, such that 



K(Px,Py,X,Y)^Yj'^"' 

n=0 

with the prescription (K in 'normal form') 
{Ko,K}^0. 



(7) 



(8) 



In these and subsequent formulas we adopt the convention of 
labeling the first term in the expansion with the index zero: in 
general, the 'zero order' terms are quadratic homogeneous poly- 
nomials and terms of order n are polynomials of degree n + 2. 
The zero order (unperturbed) Hamiltonian, 



1 



1 



Ko = Ho = -{Pi + X') + —{Pi, + Y% 



(9) 



with 'unperturbed' frequencies wi - 1,(^2 - 1 /^, plays, through 
the fundamental equation ([8]), the double role of determining the 
specific form of the transformation and of assuming the status of 
second integral of motion. 

The generating function of the transformation is a series of 
the form 



G = Gi +G2 + .. 



(10) 



and, since the procedure is based on working at each order with 
quantities determined at lower orders, the normalization algo- 
rithm proceeds by steps up to the 'truncation' order A^. At each 
step n (with 1 < « < A^), the series are 'upgraded' express- 
ing them in the new variables found with the normalizing trans- 
formation. In the Appendix A we detail the expression of the 
normal forms and the generating function we will need in the 
following. 

It is customary to refer to the series constructed in this way as 
Birkhoff normal forms. The presence of terms with small denom- 
inators in the expansion, forbids in general their convergence. It 
is therefore more eff'ective to work since the starting point with 
resonant normal forms (Sanders, Verhulst & Murdock 120071 
Gustavson 1 19661 1. which are still non-convergent, but have the 
advantage of avoiding the small divisors associated to a particu- 
lar resonance. To catch the main features of the orbital structure, 
we therefore approximate the frequencies with a rational num- 
ber plus a small 'detuning' (Contopoulos & Moutsoulas 1 1 9661 
de Zeeuw & Merritt [T983T l 



W2 



mi 

OT2 



+ d. 



(11) 



We speak of a detuned {m\:m2) resonance, with m\ + mi the 
order of the resonance. Each resonance allows us to describe 
a set of possible periodic orbits appearing in the dynamics: we 
have the 1:1 'loop', the 1:2 'banana', the 2:3 'fish' and so forth 
(Miralda-Escude & Schwarzschild |1989] l. Each of them, if sta- 
ble, is surrounded by a family of quasi-periodic orbits usually 
inheriting the same nickname. 

A very conservative strategy can be that of truncating at the 
lowest order N^ia adequate to convey some non-trivial infor- 
mation on the system. In the resonant case, it can be shown 
(Tuwankotta & Verhulst I2000I I that the lowest order to be in- 
cluded in the normal form in order to capture the main effects of 
the mi'.mi resonance with double reflection symmetries is 



A^r. 



2 X (mi + m2 — 1). 



(12) 



Using 'action-angle'-//fce variables /, defined through the 
transformation 



X = V2^icos0i, y = 
Px = y/Zhsinei, Py 



2J 2 cos 62, 
272 sin 02, 



(13) 
(14) 



the typical structure of the doubly-symmetric resonant normal 
form truncated at N^ia is (Sanders, Verhulst & Murdock |2007l 
Contopoulos . 2004 J 



K 



mi 



J I +m2J2+ Yj 'P^''\JuJ2) + 



k=2 



,7"- 7"" cos[2(m26'i - mi6»2)]. 



(15) 



where !P*** are homogeneous polynomials of degree k whose co- 
efficients may depend on 6 and the constant a„^„,A.q) is the only 
marker of the resonance. In these variables the second integral is 



£ = mi7i -H m272 



(16) 
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and the angles appear only through the resonant combination 

l/r = 11120] - m\62- (17) 

For a given resonance, the last two statements remain true for 
arbitrary A^ > Nmm- Introducing the variable conjugate to ij/. 



% - nixJi - m\J2, 



(18) 



the new Hamiltonian can be expressed in the reduced form 
Kili, i/'; £, q), that is a family of 1-dof systems parametrized by 
& and S. 

We are interested in the solution of the equations of motion. 
For a non-resonant (BirkhofF) normal form the problem is easily 
solved: the coefficient fl,„,,„, vanishes, so K lacks the term con- 
taining angles. Therefore the / are 'true' conserved actions and 
the solutions are 



X{t)-- 
where 



Y{t) = ^12J2 cos{K2t + Oo), 



(19) 



(20) 



is the frequency vector and Oq a suitable phase shift. 

In the resonant case instead, it is no more possible to write 
the solutions in closed form. It is true that the dynamics de- 
scribed by the 1-dof Hamiltonian K{'R,\l/\&,q) are always in- 
tegrable, but, in general, the solutions cannot be written in terms 
of elementary functions. However, solutions can still be worked 
out in the case of the main periodic orbits, for which J , 6 are 
again true action-angle variables. There are two kinds of peri- 
odic orbits that can be easily identified: 

1 . The normal modes for which one of the J vanishes. 

2. The periodic orbits in general position characterized by a 
^eii relation between the two angles, m29\ - m\92 = Oq. 

In both cases, it is easy to check that the solutions retain a form 
analogous to ( fT9] l with known expressions of the actions and the 
frequencies in terms of £, q and such that ki /k2 - mi /m2. 

By using the generating function ( fTOl ). the solutions in terms 
of standard 'physical' coordinates can be recovered (a part for 
possible scaling factors) inverting the canonical transformation 
defined by ( |A.3l l and jAAj . As discussed in the Appendix, the 
expansion ( fTOl i is actually composed of even-order terms only. 
Since in our applications we will treat the 1:1 and 1:2 symmetric 
resonances, we have from ( fT2b that at most A^min = 4 so that the 
transformation back to the physical coordinates expressed as a 
series of the form 

X{t) = Xi + X2 + X3 + ... (21) 

is explicitly given by 

jci = Z, (22) 

X2 = 0, (23) 

JC3 = L2(X)^{G2,X}, (24) 

X4 = 0, (25) 

JC5 = L4(X) + iLJiX) = {G4, X] + HG2, {G2, X}}. (26) 

We again remark that the vanishing of terms of even degree is 
related to the double reflection symmetry embodied in the nor- 
mal form. From a knowledge of the normalized solutions ( fT9l ), 
we can therefore construct power series approximate solutions 
of the equations of motion of the original system 

d^x 
If 



We are investigating a non-integrable system. This implies 
that any perturbation approach to cope with its dynamics is 
deemed to fail, since it produces series which do not converge in 
general. On the other hand, what the normal form provides us is 
an efficient way to construct series with an asymptotic character: 
this means that at some point we should reach an 'optimal' value 
for the expansion order A^opt (hopefully > A^min) which gives the 
best possible result (Efthymiopoulos et al. 2004). The optimal 
order depends on the size of the phase-space region we are in- 
terested in. The bigger this region is, the lower are the value of 
A^opt and the accuracy of the approximation. In galactic dynam- 
ics (contrary to what happens in celestial mechanics) it is usually 
preferable to get an overall picture of the dynamics giving up 
extreme accuracy, so that truncating at N^^^ seems a reasonable 
choice. To verify if this conjecture is actually tenable is another 
aim of the present work. 



3. Axial orbits 

In systems of the form (|6|l the orbits along the symmetry axes 
are simple periodic orbits. It can be readily verified that these or- 
bits correspond to the two normal modes for which either 7i or 
J2 vanish. If the axial orbit is stable it parents a family of 'box' 
orbits. A case that is both representative of the state of affairs 
and useful in galactic applications is that of the stability of the 
jc-axis periodic orbit (the 'major-axis orbit', if q is in the range 
(|2|i). Among possible bifurcations from it, the most prominent is 
usually that due to the 1:2 resonance between the frequency of 
oscillation along the orbit and that of a normal perturbation, pro- 
ducing the 'banana' and 'anti-banana' orbits (Miralda-Escude & 
Schwarzschild [T989b . Therefore, to get explicit solutions both 
for the major-axis orbit and the stable bananas (the 'pendulum- 
like family' in the denomination of CS90, see Section 5 below) 
we use the 1:2 symmetric normal form. 

From the expression of K reported in the Appendix, we get 
on the normal mode J2 - 0, 



3 
4' 



Ka ^2qJy-'-qj\ + -q\^-^q(q-\)\j\. 



17 



(28) 



The value of the action can be computed by using the rescaling 
dA.ll ), namely K^ - 2qE and inverting the series. E is the origi- 
nal 'physical' energy and can also be expressed by means of the 
amplitude A of the axial orbit 



£ = ilog(l+A2). 

The frequency is given by the usual differentiation 
1 dKA 



Kl 



2q dJi ' 



(29) 



(30) 



where the rescaling of the energy is taken into account in order 
to be able to use t as the physical time. Therefore, in the normal- 
ization variables we have a solution of the form ( fT9] l with F = 
and (Belmonte et al. |2007l l 



^1 



K\ 



E + -E^ 



25 . 

i 

192 



11 



!__£■ + —£"■ 

4 64 



(31) 
(32) 



-v.y. 



(27) 



Inserting this solution in the transformation formulas 

and exploiting the terms of the generating function of ( |A.16t 
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Table 1. Relative energy variations along the major-axis orbit 
with different analytic predictions. 



F A r'^"** r<^' /'' /'' 

C /\ A^p X^f Jl.(.p Ji^p 


0.1 


0.47 


0.0013 


0.0005 


0.002 


0.0001 


0.2 


0.70 


0.005 


0.005 


0.009 


0.001 


0.3 


0.91 


0.011 


0.020 


0.022 


0.003 


0.4 


1.11 


0.017 


0.058 


0.046 


0.008 


0.5 


1.31 


0.02 


0.13 


0.08 


0.02 


0.6 


1.52 


0.03 


0.25 


0.16 


0.03 


0.7 


1.75 


0.08 


0.42 


0.29 


0.05 


0.8 


1.99 


0.15 


0.70 


0.55 


0.07 


0.9 


2.25 


0.25 


0.95 


— 


0.09 


1.0 


2.53 


0.35 


— 


— 


0.12 


and (IA.171I, after a little of computer alaebra we set 


xi = Acqs Kit, 


A' _ . 


X3 = :;x(cos/fi 


t -co 


s3/(-iO, 









32 

A^ ( 59 ^ 11 , 

Xs = COS/fif + COSJ/fi? H cos5/fif 



64 \ 48 



48 



(33) 
(34) 

(35) 



This result coincides with that obtained by Scuflaire d 19951 1 
with an independent approach based on the Poincare-Lindstedt 
method and provides the explicit time evolution of an oscillation 
starting at rest from x(0) = A. 

To evaluate the quality of the approximation, a simple 
method is to follow the energy variation along the solution in 
the true potential ([TJ. We therefore compute 



1 (dxy 1 



E(^^--2\-^ --2^-^(^-<'r^ 



(36) 



and compare it with the given value of E fixed by ( |29] l for various 
amplitudes. To understand the question of the optimal order we 
can choose two different truncations of the prediction obtained 
with the normal form: 



v.(3) 






X\ + Xj,, 

xi + X2 + xs 



(37) 
(38) 



(39) 



and compute the quantity 

AE _ E(t) - E 
~ ~ E ■ 

In Fig. [T] we plot AE/E for £ = 0.1 (corresponding to an 
amplitude A - 0.47) over a half period: the curves repeat them- 
selves in the subsequent half period. The solid line is computed 
with x^j^ and the dashed line with xj^]^: the relative error in the 
energy conservation is almost three times smaller with the higher 
truncation and as low as 0.05%. However, from Fig.|2]we see that 
with E - 0.5 (A - 1.31) the situation is upset: the lower order 
truncation (which corresponds just to the first non-zero term in 
the normal form) gives an error at least five times smaller than 
that with the higher truncation. We deduce that somewhere be- 
tween the two energy levels the optimal order decreases by two 
and verify that, to get informations about an orbit 3 times bigger, 
we must accept a relative error of a few percent. In fact, in Table 
([T]i, we list the maximum absolute energy variation over a half 
period for various values of E and see that the optimal order is 




Fig. 1. Relative energy error along the major-axis orbit with two 
different truncations of the normal form at £ = 0.1. 



0.02 - 




/~\ ^\ ^^V 




^\ \ 


/ v^ \, / / 


-0.02 - 


N. \ 


u^ /^~-.~^ 


-0.04 - 




\ 1 


-0.06 - 
-0.08 - 




\ 1 
\ / 

\ / 


-0.1 - 




\ / 


-0.12 - 




\ ^-^ / 



Fig. 2. Relative energy error along the major-axis orbit with two 
different truncations of the normal form si E - 0.5. 



> 4 up to £■ = 0.2: for greater values of the energy the optimal 
order is just 2, namely with x^^ we have the best result. 

Once reached the optimal order, it can be disappointing to 
discard terms coming from a costly high-order computation. 
There are however other rules for 'summing' divergent series 
which make use of all terms (Bender & Orszag, ,1978) . like the 
construction of Pade approximant. A related approach is that of 
constructing continued fractions: successive approximants ob- 
tained by truncating the fraction at various order may give an 
improvement in the asymptotic convergence with respect to the 
original series (Khovanskii, 1963). From the normal form series 
( I37I38I ) we may compute the truncated fractions 



(3) 
.(5) 



x\ 



1 - x^jxi ' 
x\ 



1 



_i3/£l 
—J— 



(40) 
(41) 



These approximations produce rational solutions and is therefore 
natural to think to a relation with the Prendergast-Contopoulos 
approach of CS90. By using the explicit forms of ( [33H35] l we get 



.(3) 



^CF 



Acos/cif 



1 -I- jg(l -cos2a:i/) 

l+A2(|| + icos2^iO 



= Acos/fi? 



l+^'(i + l^cos2^if) 



(42) 



(43) 
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with phase shift Oq - n/l. The actions and frequencies can be 
obtained from the following procedure: starting from the normal 
form ( IA.8J4A.9] l, we determine the fixed points of the reduced 
Hamiltonian KCR, t//; S, q) with 



l/r = 01 - 6*2, 
n ^ Jx- Jl. 



(47) 
(48) 
(49) 



The fixed point corresponding to the loop is given by i/' = njl 
and 



Fig. 3. Relative energy error along the major-axis orbit with two 
different truncations of the continued fraction at £ = 0.5. 



Ji(L)i&,q) = {& + 'Rd&,q))l2, 
J2(L)i&,q) = {&-'RL(&,q))/2, 

where HiiS, q) is the solution of the algebraic equation 

dK 



dn\ 2 ^1 



(3) 

and it turns out that the expression of x^p has the same structure xhe frequency is then given by 
of the trial rational approximation used in CS90: 



A cos K\ t 
1 + Bcos2K\t 



(44) 



v.(3) 



The main difference is that in x\,p the parameters A and k\ are 
known in analytic form by ( |29] l and ( |32] i. whereas in (l44l i. A, B 
and K\ have to be numerically computed by solving a nonlinear 
algebraic system obtained by inserting the trial solution into the 
equations of motion. 

In Fig. [3] we plot the same quantities of Fig.|2]now obtained 
via the continued fraction truncations: the solid line is computed 
with x\Jp and the dashed line with xj^^. At this energy level, the 
prediction with x^' starts to outperform x),p. From Table [T] we 



CF '^'■"^'■'' "-^ v/ui.p\.^±±v-'±A±A NF' 

see that the performance of x\,p is the best when going to higher 
energies and is at least as good as that of CS90 in the same range 
of energy. 

4. Loop orbits 

As a first example of boxlet we treat the 'loop' orbits for which, 
to get explicit solutions, we can use the 1:1 symmetric nor- 
mal form. For moderate ellipticities (q > 0.7), loops ensue 
as the lowest energy bifurcation due to the 1:1 resonance be- 
tween the frequency of oscillation along the short (y-axis) peri- 
odic orbit and that of a normal perturbation (Miralda-Escude & 
Schwarzschild |1989] l. From (fT2] ). for the 1:1 resonance we have 
A^min = 2 so that a normal form truncated at K2 is already able 
to produce loops. The bifurcation curve in the (q, £')-plane starts 
from the point (1,0) (Scuflaire 1995t Belmonte et al. i2007) and 
can be expressed as the series 



E,(q) ^2(1 - q) + {I - qf --(I- q)\ 

O 



(45) 



if the normal form is truncated at progressively higher orders. 
We limit ourselves to the case q = 0.9 with transition energy 
Ec{Q.9) - 0.21 and investigate the analytic prediction of the the- 
ory by fixing the energy level at £ = 1 : with suitable rescaling, 
these are the same values of the parameters used in CS90. 

In the normalization variables we have a solution of the form 
(Owith 



(46) 



1 dK 



kl = - 



q dJ 



I(L) 



(«4;£,,)^ 



(50) 
(51) 



(52) 



(53) 



Explicit expressions of actions and frequencies are (Belmonte et 
al. l2007l l 



J\(L){&,q) 
J2(L)i&,q) 



(3-q)& + Aq{q-l) 

3q^-2q + 3 
q(3q-l)6-4q(q-l) 
3q^-2q + 3 



'-IH^ 



(54) 
(55) 
(56) 



For the solutions in the physical variables, we first work out 
the transformations (I22l - l24l i with the generating function ( lA.llI) 
obtaining 



xi - ^J2 Jul) cos KLt, 

X3 = -T y/ZJuL) cos Kit X 

OihiL) + qJl(L)) - 2{2J2(L) + qJ\(L)) COS IklI) , 

y\ - yJThiL) sin Kit, 
1 



^3 = --—ylZh(L) sin Kit X 
loq 

OihiL) + qJ\(L)) + 2(J2(L) + 2qJi(L))COS IklI) . 



(57) 

(58) 
(59) 

(60) 



We observe that, akeady in the first higher-order terms both ac- 
tions appear, to testify the strong coupling between the degrees 
of freedom. The inversion of the series 



E = -KdE,,q) = -K('Rd£>,q),'^\£',q] 
q q ^ II 



(61) 



allows us to express actions and frequencies in terms of the phys- 
ical energy. However, using the exact solutions ( 1541455] ) to evalu- 
ate ( l57l460] l would result in messy expressions hindering the pro- 
cedure of resummation with the continuous fraction. Therefore, 
in analogy with the series written for the axial orbit, a separation 
of terms of given orders is necessary and is obtained by linearly 
expanding the actions in the form 



X{t) - yjlJi COS Kit, Y{t) - V2/2 sin Kit, 



J\(L) = a(E-Ec), 
Jkl) = b + c(E-Ec). 



(62) 
(63) 
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Fig. 4. An orbit of the loop family at £ = 1.0 for q - 0.9: dots 
correspond to the numerical solution; the continuous line cor- 
responds to the prediction given by the continued fraction trun- 
cated at order 3; the dashed line that given by the normal form 
truncated at order 3 



The first expansion provides an approximate value of 7i above 
the bifurcation energy and clearly must be considered zero for 
E < E^. The second one joins in E - E^ with the corresponding 
expression on the normal mode. Inserting these into the solutions 
above and expanding in powers of E - E^ we are able to group 
terms according to their order. Clearly, this grouping does not 
affect the series themselves (namely x]^p and yj^p) but rather it 
influences the computation of the truncated fractions: we get 



(3) 
"■CF 



„(3) 



- 5 yj2a{E - £,c) cos/c^f , 



ycF = b^'^^iriKit 
where 



A\ + A2 cos iKit 
(72 H- 35/7+ 10/7 COS 2/fz,0^ 
B\ + B2Cos2/f/,f 



(64) 
(65) 



Ai = 4 (160 -I- 70/7 - (63fl -1- 70c)(£ - EJ) , (66) 

A2 = -8{20b-(9a + 20c)(E-E,)), (67) 

Bi = 18 ^(2/7(72 H-35fe)-3(21fl/7-t-35cfeH-34c)(£-£c)068) 
B2 = 36V2^(10^-3(6flH-5c)(£-£J). (69) 

For moderate values of the bifurcation energy, corresponding to 
large values of q in the range (|2]l, a simple approximation is given 
by the linear term in (|45j, E^. = 2{l - q). At this level of approx- 
imation the constants appearing in the above solutions are 




Fig. 5. The same orbit of the previous figure (dots) compared 
with the predictions truncated at order 5 (continued fraction, 
continuous line; normal form, dashed line). 




Fig. 6. Relative energy error along the loop orbit with two differ- 
ent truncations of the continued fraction. 



2(1 

2 



■q)^E^ 



(71) 
(72) 



a- -^-q. 



(70) 



and can be used to plot the orbits and compare them with numer- 
ical computations. With the choice of the parameters mentioned 
above, in Fig.|4]we compare a numerical computation of the loop 
orbit (dots) with the analytic predictions given by Xj^p (dashed 

line) and x)^^ (continuous Ime). It appears clear how the rational 
solution coming from the continued fraction truncated at order 
3 is overall quite accurate in locating both the shape and the ex- 
trema of the orbit and overtakes the prediction with the standard 
truncated series. In Fig.|5]we compare the numerical computa- 
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tion of the same orbit (dots) with the analytic predictions given 
by Xj^p (dashed line) and x^^- (continuous line). The prediction 
with the higher order truncation does not appear so much better 
as far as the shape of the loop. However, in Fig. |6] we may com- 
pare the plots of AE/E versus time for the two truncations (order 
3, continuous line; order 5, dashed line) over a half period: the 
relative eiTor is some 25% at order 5 contrary to more than 40% 



(3) 



with jc 



5. Pendulum-like (banana) orbits 

The bifurcation of the banana orbit from the major-axis occurs 
along a curve in the {q, £')-plane starting from the point (1/2, 0) 
(Scuflaire l 19951 Belmonte et al. ,2007j . It can be expressed as the 

series 

E,iq) = 8 (g -0 - y (q-if + ^iq 'if - (73) 

As before we adopt the same values of the parameters (up to a 
suitable rescaling) used in CS90: the ellipticity is ^ = 0.6 with 
transition energy Ec{0.6) - 0.76 and the energy level is fixed at 
£= 1.15. 

In the normalization variables we have a solution of the form 
(dill with 

)SKBt, Y{t) ^ ^/2J2 COS IkbI. (74) 

To find actions and frequencies, starting from the normal form 
( IA.131JAT5] i, we determine the fixed points of the reduced 
Hamiltonian KCR, i^; S, q) with 

& = Ji+ 272, 
(A = 201 - 02, 
■R = 27i - 72- 

The fixed point corresponding to the banana is given by i^ = 
and 

JuB)(&,q) = i& + 2'RB(&,q))/5, 

J2(B)iS,q) = i26-'RBiS,q))/5, 

where 'RBiS', q) is the solution of the algebraic equation 

d£. 

The frequency is then given by 
1 dK 



i'R,0;&,q) = 0. 



(80) 



Kb 



-i'RB,0;&,q). 



2q dJi(B) 

"RBiS, q) and, as a consequence 7i(B), 72(B) and kb are quite cum- 
bersome algebraic expressions involving & and q. However, sim- 
ple expressions to represent orbits in the initial physical coordi- 
nates can be obtained by replacing them with some suitable se- 
ries expansions. For, we first work out the transformations (|22| - 
with the generating function ( IA.16l4ATT7b obtaining 



yi 
yi 



727^^cos/(-Bf, 

— V27i(B) cos KBt(l(2J2{B) + ^qJ\(B)) 

-(472(B) + 6qJnB)) COS IkbI + 272(b) cos AkbI^, 
^|2J^)S\n2KBt, 

— ^J2J2(B){24■qJ^B) + 6(372(B) + 2qJnB))cos2KBt 
-8^7i(B) cos AkbI - 372(B) cos 6KBt) 



(82) 

(83) 
(84) 

(85) 



and analogous expressions for xj and 3^5. In analogy with the 
procedure followed for the loop orbit, a separation of terms of 
different low orders is useful and is obtained by linearly expand- 
ing the actions in the form 



J\(B) - a[) + a](E - Ec), 
72(B) = bi(E-E^). 



(86) 
(87) 



Inserting these into the solutions above and expanding in powers 
of E - Ec we are able to group terms according to their order. 
Here again, this grouping does not affect the series themselves 
(namely xj^^ and yj^^, with A: = 3, 5) but rather it influences the 
computation of the truncated fractions, namely the expressions 
( |40l44n i and analogous for the y coordinate. 

For moderate values of the bifurcation energy (and of orbital 
energy), corresponding to small values of q in the range (|2|, a 
simple approximation is given by the linear term in ( |73l ), E^ = 
8(^- 1/2), so that 



Jm = \E + iq-i)U+^E], 



11(B) 



-E-{q-i)\2--E 



37 



(88) 
(89) 



In this way, the expansions can be written as series of the form 






^AijCOS{2j- \)KBt + 



;=i 



(75) 




(E - he) cos KbI > Aij COS 


(76) 




j=0 


(77) 


yfr 


= Y^ A2j COs2j KBt + 


= 




5 


(78) 




(E - Ec) 2, Mj cos 2j kbI, 


(79) 




;=o 



(90) 



(91) 



so that, using ( |40l ), one can also construct jcLi. Analogously, 



(5) 



we may proceed with the higher-order truncations x)^p from 

which to obtain x^Jp. A comparison of the structure of these 
predictions with the rational solutions based on the Prendergast- 
Contopoulos approach shows that they have the same parity 
in the trigonometric parts: although in the expansions ( |90y9TT i 
many more harmonics appear, this is clearly not necessarily an 
indication of greater accuracy. 

With the choice of the parameters mentioned above, in Fig.|7] 
we compare a numerical computation of the banana orbit (dots) 
with the analytic predictions (continuous lines) given by x)^p and 

x\,p. This one, the rational solution coming from the continued 
fraction truncated at order 3, is characterized by a pair of singu- 
larities in Jf^pit) due to the presence of poles. However, the pre- 
diction is overall quite accurate in locating both the shape and 
the extrema of the orbit and overtakes the prediction with the 
standard truncated series. In Fig. |9] we plot the coiTesponding 
AE/E (continuous line) over a half period: the abrupt increase 
of the relative error is evident at the poles of the solution. 

In Fig. [8] we compare the numerical computation of the same 



,(5) 



,(5) 



orbit (dots) with the analytic predictions given by Xj^p and x 
The two predictions now almost overlap but it could be seen a 
better performance of the continued fraction truncation at the 
extrema of the orbit. In Fig.|9]we plot the coiTesponding AE/E 



Giuseppe Pucacco Dino Boccaletti and Cinzia Belmonte: Periodic orbits in the logarithmic potential 







- 








- 


1 ^ 


''^^'''''-^s^-^ J 


^ 


■^ 


^ 


^: 



Fig. 7. An orbit of the pendulum-like (banana) family at £ = 
1.15 for q - 0.6: the dots correspond to the numerical solution; 
the continuous lines correspond to the predictions truncated at 
order 3. 




Fig. 8. The same orbit of the previous figure (dots) compared 
with the predictions truncated at order 5 (continuous lines); for 
more clarity, the y-scale is expanded. 




Fig. 9. Relative energy error along the banana orbit with two dif- 
ferent truncations of the continued fraction. 



(dashed line) over a half period: the relative error is now less 
than 20% contrary to the 30% with x-'^^.. 

A comparison with the results of CS90 is possible only for 
what concerns the reconstruction of the shape and location of 
the orbit (we have used the same values of the parameters, when 
properly rescaled) and we can deduce an accuracy of our ana- 
lytic predictions at least as good as that in CS90. There is no 
information in CS90 about the ability of their solution in con- 
serving energy. 



ing line has been that of exploiting normal form expansions trun- 
cated to the first order incorporating the resonance correspond- 
ing to the given family of periodic orbits. In this way, analytic ap- 
proximate solutions can be worked out with a complete algorith- 
mic procedure. Although the attempts have always been made 
by truncating series to the first non-trivial orders, the solutions 
are definitely simple only in the case of the axial orbits (normal 
modes). For the low-order boxlets (loops and bananas), even the 
truncations at the first non-trivial order are quite cumbersome 
and would require the use of an algebraic manipulator However, 
further simplifications can be attained if the algebraic solutions 
giving actions and frequencies are expanded around energy and 
ellipticity corresponding to the bifurcation of the family. In this 
case, quite simple compact expressions of the expansions can be 
worked out, both as standard series and as continued fractions. 

A comparison with the rational (Prendergast-Contopoulos) 
approach, Contopoulos (120041) . allows us to state the follow- 
ing conclusions: the two methods are almost equivalent for what 
concerns the precision of the analytic prediction when performed 
to the same order (Contopoulos & Seimenis'1990). However, the 
normal form perturbation expansions, even if computationally 
heavy, are completely algorithmic and analytic at every stage, 
whereas the explicit evaluation of the coefficients in the rational 
expansions require the numerical solution of non-linear systems. 
We have shown that the analytic rational solutions obtained in 
this way offer a degree of reliability such that both loops and ba- 
nanas are quite well reconstructed in shape and dimension. We 
have extended the analysis in Contopoulos & Seimenis (119901 ) 
also to check the energy conservation along the boxlets: it turns 
out that, whereas for normal modes energy is conserved within 
a few percent, for loops and bananas, at this level of approxima- 
tion, it is not easy to go below 10%. 

On the theoretical side, the general usefulness of rational so- 
lutions can be explained in the light of the better convergence 
performance of truncated continued fractions. These come into 
play as a resummation method of the series expansions produced 
in the usual way in the normalization approach. The generality 
of this setting allows us to envisage analogous results in the case 
of higher-order resonances and the corresponding higher com- 
mensurable boxlets. 

In addition to the formal and algorithmic improvements, we 
remark on the relevance of this work also in relation with spe- 
cific problems of galactic dynamics. The study of orbits in non- 
axisymmetric potentials is usually performed numerically; how- 
ever, an exhaustive study with conventional integration methods 
is costly and difficult to interpret (Touma & Tremaine ll997l l. 
The availability of simple and accurate analytical recipes can 
be quite useful in several contexts in which periodic orbits and 
boxlets play an important role: we mention the study of the pa- 
rameter space of non-axisymmetric discs (Zhao, Carollo & de 
Zeeuw [r999l Zhao 1 199 9) and that of the orbit structure around 
massive black holes in galactic nuclei (Sridhar & Touma 19991 ). 
Even more promising seems to be the possibility of getting ac- 
curate solutions for periodic orbits in the triaxial case with and 
without rotation, for which the analysis is still at the level of the 
first-order averaging method applied to the 1:1:1 resonance by 
de Zeeuw (fT985l ). 

Acknowledgements. We thank G. Contopoulos for arousing our interest in this 
problem. 



6. Conclusions 

We have seen how to construct approximate solutions for the 
main periodic orbits in the cored logarithmic potential. The guid- 
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Appendix A: Resonant normal forms for the 
logarithmic potential 

In order to implement the normalization algorithm, the original 
Hamiltonian (|6]l is rescaled according to 



MtH 

CD2 



(A.1) 



so that we redefine the Hamiltonian as the series 



k=o 






(A.2) 



k=0 



where the detuning has been introducing as in ( fTTT ) and the po- 
tential is the series expansion given by (O. The normalization is 
performed with the technique of the Lie transform (Gerhard & 
Saha |199T] Yanguas|200l]l. Considering a generating function G, 
new coordinates P, X result from the canonical transformation 



(P,X) = Mc{p,x). 



(A.3) 



The Lie transform operator Mc is defined by (Boccaletti & 
Pucacco [T999] l 



k=Q 

where 

Mo = 1, Mk 



Z i^<^.'^^-^- 



(A.4) 



(A.5) 



The linear differential operator Lq is defined through the Poisson 
bracket, Lg() = {G, ■} and the functions Gj are the terms in the 
expansion of the generating function. It turns out that Go = 1 
so that, in practice, the first term in its expansion can be ignored 
as in ( [Tol l. The terms in the hew Hamiltonian are determined 
through the recursive set of linear partial differential equations 



K„='H„+Y,M„-j'Hj, n=l,2. 



(A.6) 



;=o 



'Solving' the equation at the «-th step consists of a twofold task: 
to find K„ and G„. We observe that, in view of the reflection sym- 
metries of the Hamiltonian ( |A.2| i. the chain ( |A.6l l is composed 
only of members with even index and so the normal form and the 
generating function are composed of even-index terms only. The 
unperturbed part of the Hamiltonian, 'Ko, determines the spe- 
cific form of the transformation. In fact, the new Hamiltonian K 
is said to be in normal form if, analogously to dHJ, 



{n),K]^0, 



(A.7) 



is satisfied. In the following formulas we list the normal form 
and the generating function for the expansion of the logarithmic 
potential in the cases of the 1 : 1 and 1 :2 resonances (Belmonte et 
al. 12007 J . The normal form is given in the more compact version 
given by using the action-angle-/;A;e variables /, 0: the resulting 
expressions are in agreement with the general structure of ( fTSl l. 
For the generating function it is more useful to write the explicit 
version in standard P, X variables that, although cumbersome, 
is that exploited in the transformation back to the original p, x 
variables. For the 1 : 1 resonance the terms of the normal form 
are 

Ko ^ Ji+ 72, (A.8) 

K2 = SJi - ^yf + 3 y2 + Kj^ + Kj^ ^os(20i - 202IA.9) 
o oq L 4 

4\3 16/ ' 192^2 2 
1 /39 



- (3 - ^j J\h - (1 - ^) -^lA' cos(20i - 
ie of the generating function are 



and those of the generating function 

Gi = -^P\X - —PxP\X - -|Pjf X"' - —PI PyY - 

32 X 32 ^ 32 32 ^ 

pIy PyX^Y PxXY^ PyY\ 

32q ^ 32 32 32^ 

G4 = —pIx — 

96 ^ 



202) (A. 10) 



(A. 11) 



^^PlX- 
256 ^ 



./=i 



— pI pIx -pI pIx - 

128 ^ ^ 192 ^ '^ 

13 4 7 4 19 4 19 , , 

PxPv^ + PxPiX PxXY^ PyX^Y^ + 

768 ^ 384^ ^ 768 384 

9 , , 37(7 , . 25 , . 19fl2 , 

PxPyX -PxPyX^ + PyX^Y^ —PxX^ - 

128 ^ 384 ^ 192^ 256 

P\PyY -P\PyY pI pIy + pI pIy + 

256 ^ 384 ^ 384 ^ ^ 192^ ^ ^ 

PxXY'^ + —PlPyX^Y - —PlPyX^Y 

384^ 64 ^ 64 ^ 384 

5 , , 23 

■pIx^y + — 



384^' 



— pIx^y - 

384 ^ 
23 4 19(7 4 9 . , 

PyX^Y -PyX^Y + PIXY^ - 

256 384 128 ^ 



10 
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PiXY^ - —PxPiXY"- + 



64 



Xfy^ 



31q 
384' 



192^^^'^' "384^^^'' 



pIPyY^ 



11 

64^ 
5 



PxP\XY^ + ^Px^^^^ 



3Uq 



-pIpvY^ + 



of considering the detuning term in (IA.2I ) of 2nd order in the 
perturbation: it can be shown (Pucacco et al. I2008I I that this 
choice, in principle not unique, is the 'optimal' one. Being 



X^Y^ 



31 



p3 y3 I 

2 i' present in G at order 4, they appear in K only at order 6 



^Apix' _ }IS.pix' + }hp,x' + ^ 

36 ^ 96 ^ 96 768^2 



|(3Pix + 5P.xV 



32 
64 



— (iPxP^X - P\PyY + PyX^Y + 9PxXY^) . 
64 ^ ^ 

For the 1 :2 resonance the terms of the normal form are 

Kq — Ji + 7.J2, 
K2 = 26Ji 



288^- 

2. The normalizing variables P, X have to be considered as new 

p Y^ + P^ y^nonical variables at each step of the normalization: so, for 

768q^ example, the P, X arguments of G2 when truncating GatN = 

2 are different from the arguments of G2 when truncating G 

at N - 4. A notation able to represent these features could be 

introduced but it would be heavy and we prefer to stay with 

(A.lflife standard practice of ignoring these subtleties. However, 

this observation gives reason for the apperance of the extra 

Poisson brackets in the transformations of the form ( 



%jU-j~ 



q 



J1J2, 



/5 17 \ , 29 , 

^^ = n6"i6r'^%?'^^ 

-qjUi cos(46li - 202) 



(A. 13) 
(A. 14) 



(A.15) 



and those of the generating function are 



3q 



1 



G2 = -T^Pix - -PxPix - -4Px X' + 



16 



32^ 
5?. 



P'yY 



Ga = -4:PyX 



3 
7 

"24 

13^ 
~64~ 



X^y^ 



5?: 



16 

1 



PyX'-Y - -PyXY 



24 
32^' 



PxPyY- 



-PyY' 



PiX + 



11 

144' 



pi pIx ■ 



53q 
T92' 



Pi PiX 



48 

4 4 1 4 13 4 17 

—PxPyX + PxPyX PxXY* . 

45 ^ 32^^ 90 2880 

65 .P^X^Y' _ 



-PyX^Y^ + 



PxPIx^ -PxPIx^ + 

144 ^ 192 '^ 768^ 



19^ 
~64~ 



(A. 16) 



PxX' 



1 
l44' 



P\PyY- 



19 q 

384 



P^PyY- 



2^^^ ^^^^^2^^^^'^^^ 



^^PxXY^ . l-fiPyX-Y - f^PlPyX^Y - ^^PlX^Y - 
^^PlX-Y . llpyX^Y - ^^PyX^Y . l-^PlXY^ - 



13 



7 



19 



—^pIxY^ - —PxPIxY^ + -^PxPIxY^ + —PxX^Y^ - 
192 ^ 60 ^ 64^ ^ 72 



"•"in v3v2 

192^^^ ^ " 2880"^"^ 



73 , , 25 , , 1 

PlPyY^ + ^T^PlPyY^ + :^^^-^i y 



768^ 



pIy^ + 



5q 



plx' 



I3q^- 



PiX' + 



llq 



PxX' + 



31 
768^2 



18 -^ 24 

f{3Pix^5PxX^)^ 

- [2PxP\X + 2P\PyY - 2PyX^Y + IPxXY^) . 

Concerning these formulas, two remarks are in order: 



rPyY" + 



1 



768^ 



^P\Y + 



(A. 17) 



1. The monomials associated to the detuning (namely, with S 
appearing in the coefficients) are of 2 degrees lower than that 
of the specific term of a given series. For example, in G4 (de- 
gree 6) they appear with degree 4. This is due to the choice 



